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A  SEMICLASSICAL  TRANSPORT  MODEL 
FOR  THIN  QUANTUM  BARRIERS 

SHI  JIN  AND  KYLE  A.  NOVAK  * 

Abstract.  We  present  a  one-dimensional  time-dependent  semiclassical  transport  model  for 
mixed  state  scattering  with  thin  quantum  barriers.  The  idea  is  to  solve  a  stationary  Schrodinger 
equation  in  the  thin  quantum  barrier  to  obtain  the  scattering  coefficients,  and  then  use  them  to 
supply  the  interface  condition  that  connects  the  two  classical  domains.  We  then  build  in  the  interface 
condition  to  the  numerical  flux,  in  the  spirit  of  the  Hamiltonian-preserving  scheme  introduced  by  Jin 
and  Wen  for  a  classical  barrier.  The  overall  cost  is  roughly  the  same  as  solving  a  classical  barrier. 

We  construct  a  numerical  method  based  on  this  semiclassical  approach  and  validate  the  model  using 
various  numerical  examples. 

Key  words,  multiscale  method,  semiclassical  limit,  Liouville,  von  Neumann,  quantum  barrier 

AMS  subject  classifications.  65M06,  65Z05,  81Q20,  81S30,  81T80 

1.  Introduction.  Advances  in  nanoscale  materials  fabrication  technology  have 
prompted  the  need  for  efficient  numerical  simulation  of  quantum  structures.  However, 
simulation  is  difficult  when  the  system  reacts  over  different  length  and  time  scales 
since  the  smaller  scale  usually  drives  the  accuracy  and  consistency  of  the  solution. 

Even  when  only  interested  in  the  macroscopic  behavior,  one  may  be  forced  to  resolve 
the  microscopic  dynamics.  Correspondence  principles  allow  us  to  extract  macroscopic 
behavior  from  microscopic  dynamics  in  terms  of  a  weak  limit.  When  the  scales  act  over 
several  orders  of  magnitude,  the  numerical  solution  to  the  problem  at  the  smallest  scale 
becomes  computationally  intractable.  In  these  cases,  one  often  relies  on  a  multiscale 
approach  to  provide  a  numerically  efficient  solution. 

An  example  is  the  modeling  of  electron  transport  in  nanostructures,  such  as  res¬ 
onant  tunneling  diodes,  superlattices  or  quantum  dots,  where  quantum  phenomena 
in  localized  regions  of  the  devices  cannot  be  ignored.  While  one  can  use  quantum 
mechanics  in  the  entire  region,  it  is  clearly  more  computationally  efficient  to  take  a 
multiscale  approach  using  classical  mechanics  in  the  rest  of  the  device  via  a  domain 
decomposition  technique.  Such  a  model  was  introduced  by  Ben  Abdallah,  Gamba  and 
Degond,  in  which  the  interface  conditions  connecting  the  classical  and  the  quantum 
regions  were  introduced  to  couple  two  classical  regions  with  a  quantum  region  [6,  7,  8], 

This  work  is  an  extension  of  the  Hamiltonian-preserving  finite-volume  method 
introduced  by  Jin  and  Wen  [17,  16]  for  solving  the  multi-dimensional  classical  Liou¬ 
ville  equation  with  a  discontinuous  (but  classical)  potential.  The  idea  there  was  to 
build  the  interface  condition,  such  as  used  in  [7],  that  properly  incorporates  partial 
transmission  and  reflection  information  at  the  barrier  into  the  numerical  flux.  This 
produces  a  scheme  that  connects  momenta  (velocities)  on  both  sides  of  the  barrier 
via  the  Hamiltonian  preservation  principle.  Such  a  method  is  stable  in  both  ll  and 
l°°  norms  under  a  hyperbolic  stability  condition  and  captures  sharply  the  weak  semi¬ 
classical  limit  of  the  linear  Schrodinger  equation  or  geometrical  optics  through  the 
barrier  or  interface. 


'Department  of  Mathematics,  University  of  Wisconsin,  480  Lincoln  Drive,  Madison,  Wisconsin 
53706-1338,  USA  (jineraath.wisc.edu,  novakamath.wisc.edu)  The  views  expressed  in  this  article 
are  those  of  the  author  and  do  not  reflect  the  official  policy  or  position  of  the  United  States  Air 
Force,  Department  of  Defense,  or  the  U.S.  Government. 
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The  quantum  barrier  that  separates  the  two  classical  regions  differs  from  a  clas¬ 
sical  barrier  in  that  a  quantum  wave  can  tunnel  through  a  barrier,  be  partially  trans¬ 
mitted  and  reflected  by  a  barrier,  and  resonate  inside  a  barrier.  Our  idea  is  to  solve 
the  Schrodinger  equation  (either  exactly  if  possible,  or  numerically  via  a  transfer 
matrix  method  [1,  18,  12])  inside  the  quantum  barrier  in  order  to  generate  the  trans¬ 
mission  and  reflection  coefficients,  and  then  use  that  information  in  the  interface 
condition  to  solve  the  classical  Liouville  equation  through  the  barrier,  in  the  spirit 
of  the  Hamiltonian-preserving  method  of  Jin  and  Wen.  When  the  quantum  barrier 
is  thin  (on  the  order  of  a  de  Broglie  wavelength),  solving  the  stationary  Schrodinger 
equation  suffices.  Thus,  our  first  step  is  merely  preprocessing.  Once  the  transmission 
and  reflection  coefficients  are  generated,  the  time  marching  is  based  on  classical  me¬ 
chanics.  Hence,  our  approach,  which  efficiently  handles  a  thin  quantum  barrier,  has 
a  computational  cost  similar  to  a  classical  simulation  in  the  entire  device. 

In  §2  we  review  the  correspondence  between  the  classical  and  quantum  mechan¬ 
ics.  We  then  describe  the  quantum  scattering  at  barriers  in  §3.  We  propose  the 
semiclassical  model  and  its  numerical  discretization  §4.  In  §5  we  present  four  numer¬ 
ical  examples  to  verify  and  validate  the  semiclassical  model  the  numerical  method. 
Our  numerical  results  indicate  that  the  model  captures  correctly  the  solution  of  the 
Schrodinger  equation  in  the  entire  domain  in  the  limit  of  vanishing  Planck  constant. 

2.  Correspondence  between  classical  and  quantum  mechanics. 


2.1.  From  classical  to  quantum  mechanics.  A  typical  problem  under  con¬ 
sideration  is  particle  flow  through  a  solid-state  device.  If  the  potential  is  sufficiently 
smooth  we  may  describe  non-interacting  particle  dynamics  in  phase  space  classically 
as  a  Hamiltonian  system 

(2.1)  ^  =  Z-  =  VpH(x,p),  %  =  -VZV  =  -VxH{x,P) 

at  TTi  at 

where  x(t )  e  is  the  particle  position,  p(t)  €  Rd  is  the  momentum,  m  is  the  effective 
mass  and  V(x)  is  a  time-independent  potential.  The  Hamiltonian  function  H(x,p) 
represents  the  total  energy  of  the  system 

(2.2)  H(x,p)  =  ^  +  V(x)=E. 


One  may  introduce  a  probability  distribution  of  particles  f(x,p,  t )  in  phase  space. 
By  requiring  that  the  probability  be  conserved  along  the  particle  trajectories  one  has 


lf  = 

dtJ 


d  dx 

Wr+lt 


■f  +  !'V,/  =  o 


and  with  the  help  of  equation  (2.1),  one  gets  the  classical  Liouville  equation 


(2.3)  ^f  =  {HJ}  =  VpfVxH-Vxf-VpH 
where  {  ■ ,  •  }  is  the  Poisson  bracket.  Alternatively, 

(2.4)  ^/+p-V*/-V*F(x)-Vp/  =  0. 

By  considering  the  zeroth-order  moment  of  f(x,p,t ),  one  obtains  the  probability 
position  density  in  physical  space 

p(x,t)=  f  (x,  p,  t)  dp, 

JRd 
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which  serves  as  a  primary  observable  for  the  comparison  of  the  model. 

When  the  potential  fluctuates  rapidly  over  a  short  distance  or  the  particles  im¬ 
pinge  on  a  sharp  jump  in  potential,  the  classical  description  fails  to  capture  the 
quantum  wave-like  nature  of  the  particle  and  the  Liouville  description  produces  an 
incorrect  solution.  In  particular,  the  classical  Liouville  equation  does  not  model  bar¬ 
rier  tunneling,  probabilistic  partial  reflection  and  transmission,  or  resonance  which 
are  crucial  to  the  behavior  of  many  modern  electronic  devices. 

By  considering  Dirac  quantization,  one  has  the  formal  correspondence  between 
the  classical  quantities  and  the  quantum  operators 

(2.5)  x  — >  x,  p  ►  —  fftV,  and  E  — > 

where  h  is  Planck’s  constant.  Using  this  quantization,  one  obtains  the  Schrodinger 
equation  from  the  classical  Hamiltonian  (2.2) 

(2.6)  ih^  ~  ^  =  (~^m  A  +  V ^ 

which  describes  the  time  evolution  of  the  probability  amplitude  ip(x,t',x,p)  initially 
centered  at  x  with  an  initial  energy  state  E  =  H(x,p).  The  square  of  the  magnitude 
of  the  probability  amplitude  p(x,t)  —  \ip(x,t)\2  gives  the  position  density  in  physical 
space. 

Instead  of  considering  a  pure  state  system,  one  may  also  consider  a  mixed  state 
system  for  which  the  initial  state  H(x,p)  of  the  particle  is  given  in  terms  of  a  macro¬ 
scopic  statistical  distribution  f(x,p).  Define  the  density  matrix  as 

(2.7)  p(x,x',t)  =  /  f(x,p)'ip(x,t-,x,p)4>(x',t\x,p)dxdp. 

J R'i  J&* 

The  time  evolution  of  the  density  matrix  is  found  by  taking  the  partial  derivative  of 
equation  (2.7)  with  respect  to  t.  By  using  the  Schrodinger  equation  (2.6)  and  the 
hermicity  of  Hamiltonian  operator  H ,  one  obtains  the  von  Neumann  equation 

(2.8)  ih^p{x,x',t)  =  “  Ax<]  -\-V{x)  -V(x')Sj  p{x,x',t). 

The  von  Neumann  representation  may  be  thought  of  as  the  fundamental  descrip¬ 
tion  of  quantum  mechanics  [9],  By  taking  f(x,p)  =  5(x  —  xq )6(p  —  p o)  in  (2.7),  the 
density  matrix  reduces  to  p(x,x',t)  —  ^(x,t\xo,po)il}{x' ,t\XQ,po)  and  the  physical 
observables  of  the  mixed  state  von  Neumann  equation  correspond  to  those  of  the 
pure  state  Schrodinger  equation.  In  this  manner,  the  Schrodinger  equation  is  simply 
a  limiting  case  of  the  von  Neumann  equation.  By  taking  the  diagonal  of  the  density 
matrix,  one  gets  the  position  density  in  physical  space 

p(x,x,i;)=  /  f{x,p)\ip(x,t-,x,p)\2  dxdp. 

Jr*  Jud 

2.2.  Semiclassical  limit:  quantum  to  classical.  Consider  a  characteristic 
length  and  time  scale  LSx  and  Ldt  where  8x  is  the  natural  length  scale  such  as  a 
de  Broglie  wavelength  5x  =  h/p  for  some  momentum  p.  By  rescaling  x,x'  and  t 

x  t~*  x/L5x ,  x'  i — *  x1  /L5x,  t  >— >  t/L5t 
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in  the  von  Neumann  equation  we  have 

(2.9)  ie^p(x,x',t)  =  ~  A*'l  +  v(x)  ~  p{x,x') 

where  the  dimensionless  scaled  Planck  constant  e  =  [mL{8x)2 / <5i]— 1  h  and  the  effective 
mass  m  has  been  nondimensionalized.  Solving  the  Schrodinger  and  von  Neumann 
equations  numerically  presents  several  difficulties.  The  de  Broglie  wavelength  must  be 
resolved  numerically  to  ensure  correct  physical  observables  of  the  solution.  Typically, 
this  requires  that  the  mesh  size  Ax  =  O(e)  or  even  o(e)  with  a  similar  constraint  on 
the  time  discretization  At  [4,  24].  When  e  is  small,  computation  is  expensive  since  we 
need  to  use  0(Nd+1)  operations  to  compute  the  Schrodinger  solution  and  0(N2d+l) 
operations  to  compute  the  von  Neumann  solution  where  N  =  0(e~1)  is  the  number  of 
grid  points  in  each  space  dimension.  Because  of  such  reasons,  semiclassical  methods 
are  important  for  the  solutions  when  e  1. 

A  typical  path  to  the  derivation  of  semiclassical  limit  is  through  the  WKB  ap¬ 
proximation.  However,  the  WKB  approximation  to  the  Schrodinger  equation  fails  to 
capture  multiphase  information  beyond  caustics  [14,  29],  An  alternative  method  is  to 
use  the  Wigner  transform,  the  Fourier  transform  of  the  density  matrix, 

(2.10)  W(x,p,t)  =  f  p(x  +  \ey,  x  -  \ey,  t)e~lpv  dy. 

(27T)a  JRi 

By  applying  the  transform  to  the  von  Neumann  equation  one  has  the  Wigner  equa¬ 
tion  [31] 


%-w  +  £  •  VXW  -  ©£W  =  0 
at  m 

where  the  nonlocal  term 

&EW(x,p,t)  =  ^  /  -  [V{x+\ey)-V{x-\ey)\W{x,y,t)e-ilp-y  dy. 

with 

W{x,y,t)=  [  W(x,p,t)eip  y  dp 

J  Kd 

being  the  inverse  Fourier  transform  of  W(x,p,  t).  When  the  potential  V(x)  is  suf¬ 
ficiently  smooth,  one  recovers  the  classical  Liouville  equation  in  the  limit  as  e  — > 
0  [10,  22] 

(2.11)  f£  +  --V,/-V,y-Vp/  =  0. 

at  m 

However,  the  classical  limit  is  not  valid  at  the  discontinuties  of  the  potential  [3, 
27,  28],  where  the  potential  behaves  as  a  quantum  scatterer.  In  the  case  of  a  quantum 
barrier,  we  may  consider  a  multiscale  domain  decomposition  approach  for  a  solu¬ 
tion  [7].  In  the  next  section,  we  present  a  semiclassical  model  of  a  thin  quantum 
barrier  with  the  mixed-state  dynamics. 
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3.  Particle  behavior  at  a  quantum  barrier.  To  model  quantum  dynamics, 
we  consider  a  top-down  multiscale  approach  by  considering  the  quantum  effects  as 
local  corrections  to  the  global  classical  particle  dynamics.  In  order  to  isolate  and 
simplify  the  problem,  we  make  the  following  assumptions/limitations: 

1.  The  dynamics  are  restricted  to  one  dimension. 

2.  The  effective  width  of  a  barrier  is  O(e).  On  the  classical  scale,  this  means 

that  we  may  approximate  it  as  having  zero  width;  on  the  quantum  scale,  this  means 

that  we  may  typify  it  as  a  single  scattering  center  and  we  may  neglect  particle  dwell 

time  in  the  quantum  region  in  the  semiclassical  limit. 

3.  The  distance  between  neighboring  barriers  is  0(1)  and  hence  each  barrier 
may  be  considered  independently. 

4.  The  change  in  the  potential  VxV(x)  is  0(1)  except  at  quantum  barriers. 

5.  The  coherence  time  is  sufficiently  short  and  therefore  we  may  neglect  inter¬ 
ference  away  from  the  barrier. 

Naturally,  one  would  like  to  be  able  to  treat  a  wider  class  of  problems  including 
periodic  crystalline  domains  and  mesoscopic  barriers  for  which  s  is  nonvanishing.  We 
will  examine  corrections  and  extensions  to  these  simplifications  in  subsequent  papers. 
We  begin  with  the  Hamiltonian  system  discussed  in  §2 

d  d 

—x  =  X7pH(x,p),  —  p=  -VxH(x,p). 

Let  a  bicharacteristic  of  the  function  H(x,p )  be  the  integral  curve  ip(t)  =  ( x(t),p(t )). 
Note  that  <p(t)  may  not  be  defined  for  all  time  f  sK.  When  H(ip(t))  is  differentiable, 

(3-1}  jtH(<p(t))  =  ±x.VxH  +  jtp-VpH  =  0 

from  which  it  follows  that  the  Hamiltonian  is  constant  along  any  bicharacteristic  <p(t), 
i.e., 

(3.2)  H(tp(t))  =  const. 

Condition  (3.1)  may  be  interpreted  as  the  strong  form  of  the  conservation  of  energy, 
while  condition  (3.2)  may  be  interpreted  as  the  weak  form.  If  the  potential  V(x )  is 
discontinuous  or  not  defined  in  some  region  Q  €  Rd,  the  Liouville  equation  fails  to 
have  a  global  solution  since  VxV(Q)  is  undefined. 

The  key  idea  behind  Hamiltonian  preserving  schemes  [17,  16]  is  to  (a)  solve 
the  Liouville  equation  locally;  (b)  use  the  weak  form  of  the  conservation  of  en¬ 
ergy  to  connect  the  local  solutions  together;  and  (c)  incorporate  a  physically  rel¬ 
evant  interface  condition  to  choose  the  correct  solution.  Let  £  be  the  locally  de¬ 
fined  set  of  bicharacteristics  of  the  function  H(x,p).  By  requiring  the  Hamiltonian 
to  be  constant  along  trajectories,  we  create  an  equivalence  class  of  bicharacteristics 
M  =  {  P*  €  £  I  H(P*)  =  H(p)  }. 

Generating  a  global  bicharacteristic  is  a  matter  of  connecting  equivalent  bichar¬ 
acteristics  at  the  barriers.  If  we  consider  the  incident  and  scattered  trajectory  limits 
(x(t~),p(t~))  and  (x(t+),p(t+))  on  a  barrier  in  one-dimensional  phase  space,  then 
from  equation  (2.2)  the  scattered  momenta  are 

(3.3a)  p(t+)  =  ~p{t~) 

for  reflection  and 

(3.3b)  p(t+)  =  sign[p(f_)]\/|p(t-)|2  +  2m{V{x(t~)  -  V(x (t+))] 
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for  transmission.  Unless  \p(t~)\2  <  2m[V(x(t+)  —  V(z(f~))],  for  which  the  trans¬ 
mitted  momentum  is  imaginary,  the  conservation  of  energy  does  not  tell  us  which 
of  these  two  bicharacteristics  a  particle  should  physically  follow.  In  order  to  resolve 
the  nonuniqueness,  we  require  an  additional  interface  condition  which  we  derive  from 
the  Schrodinger  solution  across  the  interface.  By  interpreting  a  wave  function  as  a 
statistical  ensemble  of  a  large  number  of  particles  [26],  we  have  the  interface  condition 

(3.4)  f(x(t+),p(t+))  =  R(pr{t~))f(x(t+),pr(r))+T(pt(t-))f(x(t-)!Pt(t-)) 

where  T(p )  denotes  the  probability  of  an  incident  particle  being  transmitted  across 
some  region,  R(p )  denotes  the  probability  of  an  incident  particle  being  reflected,  and 
the  incident  momenta 


Pr(t  )  =  — p(£+)  and 

Pt(t~)  =  sign[p(f+)]-v/|p(t+)|2  +  2 m[V(x{t+)  -  V(x{t~))} 


come  from  equations  (3.3b)  and  (3.3a). 

We  assume  that  the  probability  of  a  particle  being  absorbed  by  the  barrier  is  zero 
and  hence  T(p)  +  R{p)  =  1.  By  defining 


W)) 


1  if  \p(t  )|2  >  2 m[V(x(t+))  -  V(x(t  ))]  and 
0  otherwise, 


i.e.,  total  transmission/reflection,  condition  (3.4)  reduces  the  classical  Liouville  so¬ 
lution  for  which  bicharacteristics  are  uniquely  determined  for  each  ( x,p ).  When 
T{p)  €  (0, 1),  i.e.,  partial  transmission/reflection,  the  bicharacteristics  are  no  longer 
unique  and  instead  we  consider  multiple  bicharacteristic  solutions. 

Every  interaction  with  a  barrier  potentially  introduces  a  reflected  and  transmitted 
solution  resulting  in  an  additional  bicharacteristic.  We  may  enumerate  the  solutions 
and  define  a  bicharacteristic  solution  to  the  Liouville  equation  as 

fk(x,p,t)  =  J  f  {x,p)pk(x,p,t\  x,p)  dx  dp 

where 

pk{x,p,t\x,p)  =  6{x{t)  —  x)8{p{t)  -p) 

is  the  fcth  global  bicharacteristic  for  H{x,p).  By  linearity  of  the  Liouville  equation  we 
may  consider  the  general  solution  as  the  superposition  of  the  bicharacteristic  solutions 

(3.5)  f{x,p,t)  =  ^2ck(H(x,p))fk{x,p,t). 

k 

where  ck(H(x,p))  is  product  of  reflection  and  transmission  probabilities  along  the  fcth 
bicharacteristic. 

Except  for  simple  solutions  such  as  the  global-in-time  solution  for  a  piecewise- 
constant  potential  or  the  local-in-time  solution  for  a  piecewise-quadratic  potential,  an 
exact  solution  cannot  be  explicitly  given.  Even  for  a  simple  discontinuous  oscillator 
the  number  of  bicharacteristics  that  need  to  be  tracked  becomes  cumbersome  in  a 
short  time  interval.  See  Fig.  3.1.  By  solving  the  model  numerically,  we  mitigate  these 
difficulties. 
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Fig.  3.1.  Particle  position  as  a  function  of  time  for  potential  V(x)  —  2\x\  —  H(x)  where  II (x) 
is  the  Heaviside  step  function.  Particle  has  initial  conditions  S(x  +  l)5(p). 


4.  A  semiclassical  approach  and  numerical  discretization. 

4.1.  A  semiclassical  approach.  When  the  quantum  barrier  is  sufficiently  nar¬ 
row,  the  barrier  may  be  modeled  using  the  time-independent  Schrodinger  equation. 
We  may  then  derive  the  transmission/reflection  probabilities  for  the  interface  con¬ 
dition  (3.4)  by  considering  the  current  density.  The  interface  condition  is  used  to 
connect  two  classical  domains  modeled  by  the  classical  Liouville  equation  (2.11). 

We  consider  an  algorithm  consisting  of  an  initialization  routine  and  a  Liouville 
solver: 

1.  During  initialization,  we  determine  the  stationary  states  at  the  barrier  by 
solving  the  time-independent  Schrodinger  equation.  The  solutions  may  be  found  by 
considering  the  barrier  as  an  open  quantum  system  [2]  outside  of  which  the  potential 
is  constant.  Typically,  this  may  be  done  by  using  a  quantum  transmitting  boundary 
method  [20],  a  spectral  projection  method  [23],  or  a  transfer  matrix  method  [1, 18, 12]. 
With  this  solution,  we  compute  the  scattering  information,  namely  the  transmission 
and  reflection  coefficients. 

2.  Following  initialization,  we  solve  the  Liouville  equation  using  a  finite  volume 
method.  As  done  in  [15]  the  interface  condition  (3.4)  is  built  into  the  numerical 
flux  in  a  framework  called  the  Hamilton  preserving  scheme.  This  yields  a  numerical 
scheme  for  which  the  stability  condition — the  CFL  condition — is  hyperbolic,  namely 
At  =  0( Ax,  A p)  with  l°°  and  ll  stability.  See  [15], 

This  approach  aims  at  capturing  the  weak  limit  of  the  Schrodinger  and  von  Neu¬ 
mann  equations  as  e  — >  0,  without  solving  the  Schrodinger  or  von  Neumann  equations 
over  the  entire  domain,  but  rather  just  at  the  quantum  barrier  and  only  in  the  initial¬ 
ization  step.  We  now  discuss  the  initialization  routine  and  the  finite  volume  routine 
in  detail. 

4.2.  Routine  initialization.  We  use  the  transfer  matrix  method  because  it  is 
robust  over  a  wide  range  of  momenta.  On  the  quantum  scale  we  decompose  a  one¬ 
dimensional  barrier  into  a  sequence  of  step  potentials  over  which  we  solve  the  time- 
independent  Schrodinger  equation  exactly.  Take  a  quantum  barrier  in  the  bounded 
region  Q  =  [aq ,  x 2]  and  take  the  potential  to  be  constant  outside  this  barrier — V (x)  = 
Vi  in  C\  =  (—00, aii)  and  V{x)  =  V%  in  C\  =  (x2,oo).  For  a  state  E  =  p2 /2m  the 
time-independent  Schrodinger  equation 

—£ip"(x)  +  2mV  (x)tp(x)  =  p2tp(x) 
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Fig.  4.1.  Approximation  of  a  potential  barrier  by  a  series  of  step  potentials.  The  effective 
transfer  matrix  M  =  Mn  *  •  •  M2  Mi  where  M  j  is  the  transfer  matrix  for  a  step  potential  at  xj . 


has  the  solution 

’  aieiKl(x-Xl)/e  +  bie-iKl(x-Xl)/et  X  £  Cl, 

(4.1)  ip(x)  =  <  ipQ,  x  £  Q 

a2e^2(x-x2)/e  +  b2e-iK2(x-x,)/e t  xec2 

where  kXi2  =  yjp1  —  2m V)  2  and  the  coefficients  ay,  a2,  bi  and  b2  are  uniquely  de¬ 
termined  by  the  boundary  conditions  at  xy  and  x2.  By  requiring  that  the  solution 
ip(x)  and  its  derivative  be  continuous,  ipQ  is  uniquely  determined  by  the  values  ay 
and  by  using  the  boundary  conditions  iPq{x{)  and  iP'q(x  1).  In  turn,  the  values  a2  and 
b2  are  uniquely  determined  by  the  values  ip q(x2 )  and  iP'q{x 2).  Since  the  Schrodinger 
equation  is  linear,  a2  and  b2  may  be  expressed  as  linear  functions  of  ay  and  by.  Hence, 
for  each  momentum  p  we  may  relate  the  solution  in  C2  with  the  solution  Cy  in  terms 
of  the  transfer  matrix  M 


ron  mJ2 
m2 1  m2  2 


)(::) 


An  arbitrary  quantum  barrier  may  be  discretized  and  approximated  by  a  series 
of  step  potentials,  for  each  of  which  a  transfer  matrix  may  be  computed  analytically. 
Specifically,  the  transfer  matrix  may  be  approximated  as  M  =  Mn  ■  ■  ■  M2Mi  with 
M j  —  D*|2  P j  D V2  where 


P  __  1  A  +  Kj/lij+y 

3  2  V1  — 


1  Kj  j  Kj  4- 1 
1  -j-  Kj  /  Kjj.y 


is  the  transfer  matrix  associated  with  a  potential  jump  V (x^ )  —  V(xj  )  and 

,  .  n  _  fexp(iAxKj/e)  0  \ 

uJ~y  0  expi-iAxKj/e) J  ' 

is  the  transfer  matrix  associated  with  the  displacement  Ax  =  Xj  —  Xj^y. 

One  may  also  express  the  solutions  in  Cy  and  C2  in  terms  of  a  scattering  matrix 
S  which  relates  the  incident  and  scattered  waves 


(ry  t2 
Ji  r2 


— m2i/m22 
A /m22 


1  /m22 
m12/m22 


where  A  =  det  M  =  m22myy  —  mi2m21.  By  considering  the  time  evolution  of  the  posi¬ 
tion  density  p(x,  t )  =  \ip(x,  t)  |2  in  the  Schrodinger  equation,  one  derives  the  continuity 
equation 


-  0  +  V  •  J  —  0 


8 


where  the  current-density  is  defined  as  J(x)  =  em  1Im  (i/>Vi/>)-  From  equations  (4.1), 
one  has  that 

(4  cn  J(r)=  f«i(l«il2-|fcil2)M  *£Ci 

1  J  [)  {^(\a2\2-\b2\2)/m,  x£C2 

where  m  is  the  effective  particle  mass.  The  positive- valued  terms  of  the  J(x )  express 
the  flux  of  right-traveling  waves  and  the  negative-valued  terms  express  the  flux  of 
left-traveling  waves.  In  particular,  for  a  wave  incident  on  the  barrier  from  the  left 
(b2  =  0),  we  have  a2  =  t\a\  and  fq  =  r\d\.  It  follows  that  the  reflection  coefficient  R\, 
the  ratio  of  the  reflected  to  incident  current  densities,  and  the  transmission  coefficient 
Tj,  the  ratio  of  the  transmitted  to  incident  current  densities,  are 


Ri  =  |ni2 


Ti  =  (/c2/«i)|ti  |2 


Similarly,  for  a  wave  incident  from  the  right 

(4.8)  #2  =  H2  and  T2  =  {k  i/K2)|f2|2- 

The  transmission  and  reflection  coefficients  are  uniquely  determined  along  a  bichar¬ 
acteristic.  It  is  clear  by  time-reversibility  that  the  transmission  coefficient  along  any 
bicharacteristic  is  independent  of  direction 


Ti(p)  =  T2  (->/pa  +  2m(V2-V1))  . 


4.3.  A  Liouville  solver.  Without  loss  of  generality,  we  shall  take  the  mass 
m  —  1  in  which  case  we  equate  the  velocity  with  the  momentum  p.  To  solve 
the  semiclassical  Liouville  equation  (2.11),  we  use  a  Hamiltonian-preserving  finite- 
volume  method  [17].  We  consider  a  uniform  mesh  in  phase  space  with  grid  points  at 
(a'i+i/2 1  Pj+1/2)  and  denote  grid  spacing  Ax  =  a;i+1/2  -  Xj_1/2  and  A p  =  pj+ 1/2  - 
Pj-1/2  i>3  €  Z.  Let  the  cell  centers  be  Xi  =  \{xi+i/2  +  ah-1/2)  and  pj  = 

\{Pj+ 1/2  T Pj — 1 72 ) -  For  convenience  of  notation,  we  shall  take  p0  =  0  and  p_j  =  —  pj. 
We  shall  consider  the  quantum  barrier  to  be  located  at  a  cell  interface  xz+ 1/2  for 
some  integer(s)  Z. 

Define  the  cell  average  over  the  cell  C,j  =  [aq_i/2,  aq+i/2)  x  [Pj-1/2, Pj+1/2)  as 


fZ  =  A^IJc.f{-X’P’tn)dxdp- 


The  finite- volume  discretization  of  the  one-dimensional  Liouville  equation  (2.11)  is 
(4.10)  /£.+1  =  /£  -  At  \pjdxj§  -  dxVidpfjj] . 

where  the  discrete  operators  dxftj,  dpfij  and  dxVi  are 

dxfij  =  (/ i+l/2,j  ~  ft-l/2,j)l ^X-> 

dpfij  =  (/i.j+1/2  -  /i,j-i/2)/Ap,  and 
Wi  =  (T^+l/2  -  K-l/2 )/&* 
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with 


f ± 

•'i+1/2, j 
fi,j+ 1/2 


l  rPj+i/2 

hm  /  f(x>P)dP> 

Jpj-x/2 


„ _  ± 

X  Xi+l/2 


_1_ 

Arc 


rx  i+1/2 

/  f(x,Pj+i/2)dx,  and 


V  i+1/2 


Upwinding  is  used  to  approximate  the  fluxes  and  /ij+i/2-  If  the  potential 

V(a;)  is  continuous  at  some  point  £j+i/2,  then  p(£+)  =  p(t_)  and  hence  = 

P+1/2  j  which  reduces  the  discretized  Liouville  equation  (4.10)  to  the  usual  upwind 
finite  volume  scheme.  At  the  barrier  xz+1/2,  special  consideration  must  be  taken. 

From  conservation  of  the  Hamiltonian  (3.3)  we  have  that  the  incident  velocity  qj 
(upwind  of  the  barrier)  for  a  particle  transmitted  with  velocity  pj  is 


9±b'  —  |;|  ^  ^(Vz+i/2  Vz+1/2 )■ 


Similarly,  the  transmitted  velocity  (downwind  of  the  barrier)  for  a  particle  incident 
with  velocity  pj  is  —q~j.  The  incident  velocity  for  a  particle  reflected  with  velocity 
Pj  is  simply  —pj.  Note  that,  whereas  —p~j  =  Pj,  in  general  —q-j  ^  qj-  Further  note 
that  by  time  reversibility  T(q~j)  =  T(j>j)  and  R(q-j)  =  R(pj)- 

The  left  and  right  limits  of  the  probability  distribution  /  in  the  cells  immediately 
downwind  of  the  quantum  barrier  are  determined  by  the  interface  condition  (3.4) 

fz+i/2,j  =  ■^(9j)/z+i/2,-j  +  T(qj)f{xz+l/2,  qj)  for  j  >  0 

f Z+1/2J  ~  R(lj)fz+ 1/2, -j  +  ^(<b')/(a;z-fi/2’ 9/)  for  j  <  0. 

The  values  for  /(zf +1/2,  Qj)  are  approximated  in  a  manner  similar  to  Scheme  II  of  [17], 
Consider  the  flux  incident  from  the  left  (qj  >  0) — the  same  treatment  applies  to  flux 
incident  from  the  right.  We  define  f(xz+i/2’<lj')  as  the  cell  average 

1  ri:+ 1/2 

(4-11)  f(xZ+l/2  'qA  =  ^pjq  Pf(xZ+l/2>P)dP 

where 

<7j±i/2  =  \]p2j±  1/2  +  2(^z+i/2  -  Vz+i/ 2)- 

The  integral  is  approximated  by  a  composite  mid-point  rule.  Since  the  limits  of  the 
integral  are  not  generally  gridpoints  in  the  p-direction,  some  care  must  be  taken.  If 
Pfc-1/2  <  Qj- 1/2  <  Qj+i/2  <  Pk+ 1/2  some  k,  then  we  take 


f(xz+\/2'Qi)  ~  fz+i/2,k  +  1jap(f Z+\/2,k) 
where  the  slope  <xp(-)  in  the  p-direction  is  calculated  using  the  Van  Leer  limiter 


(4.12) 


Ppifij)  ~ 


fij  ~  fi, 
Ap 


;l\  (h  ( IhOl.  fij 
)  \  fij  ~ 
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with  <f>{6)  =  (0  +  |0|)/(1  +  |0|)  [21],  Otherwise  pk-i/2  <  Qj- 1/2  <  <  Qj+ 1/2  < 

Pk+s+ 1/2  for  some  k  and  s,  and  we  take 


0 


+ 


(4.13)  f{xz+i/ 2’^i)  — 

p  {(Pfc+1/2  -  Qj-1/2)  [pfc/z+l/2,fc  +  5(P*+l/2  +  qj-l/2)vp(pkf z+l/2,k) 

Pk+lApf z+l/2,k+l  ^  h  Pk+s-l^pf z+l/2,k+s-l  + 

(<7j+l/2  ~  Pfc+s-l/2)  ^Pk+sf  z+l/2,k+s  2^Pk+s-l/2  +  Qj+l/2)a’p(Pk+sfz+l/2,k+s^\\' 

For  a  second-order  accurate  method  we  use  a  slope-limited  piecewise-linear  inter- 
polant  to  approximate  the  right  and  left  density  limits 


(4.14) 


■Ci^f  1/2 j  ~  fij  T  2(1  Aj)Ascrx(/ij) 


where  A j  =  |uj|At/Ax  and  the  slope  ax(-)  in  the  x-direction  is  calculated  using  the 
Van  Leer  limiter 


(4.15) 


&x{fij)  — 


fij  — 

Ax 


<t> 


fi+l,j  fij 
fij  ~  fi-lj 


Since  the  slope  crx( ■)  is  a  function  of  fi-ij,  fij  and  fi+ij  and  the  density  /  is 
not  necessarily  continuous  across  the  barrier  in  the  x-direction,  we  can  not  directly 
use  (4.14)  and  (4.15)  to  calculate  the  density  limits  at  the  barrier  interface.  Rather, 
we  first  need  to  construct  the  ghost  densities  /£  and  /|+1  across  the  barrier  using 
the  scattered  densities  at  xz  and  xz+\  based  on  conservation  of  mass.  Specifically, 
downwind  of  the  barrier 


fz+1/2  ~~  fz+i/2^Z'  fz+ii  fz+2)  a,nd  fz+i/2  ~  fz+i/2^z~^'  fz+i) 

with  ghost  densities  /£  and  fz+i  located  upwind  of  the  barrier;  and  upwind  of  the 
barrier 


fz-1/2  —  fz-\/2^z~^'  fZ’  fz+i)  ant^  f z+3/2  —  ^z+z/2^fz>  fz+i>  fz+2) 

with  ghost  densities  /£  and  fz+ 1  located  downwind  of  the  barrier. 

Construction  of  the  ghost  densities  is  analogous  to  using  ghost  cells  to  enforce 
semipermeable  inflow  and  outflow  reflecting  boundary  conditions.  To  calculate  the 
ghost  densities  upwind  of  the  barrier  we  use  the  interface  condition  (3.4)  to  mix  to¬ 
gether  the  densities  upwind  of  the  barrier  that  will  subsequently  be  combined  through 
transmission  and  reflection.  In  this  case, 


fzj  =  R{Qj)fz+i,-j  +  T{qj)f(xz,  qj)  for  j  >  0, 

fz+ij  =  R{lj)fz-j  +  T(qj)f{xz+I,qj)  for  j  <  0. 

To  calculate  the  ghost  densities  downwind  of  the  barrier  we  unmix  the  densities  down¬ 
wind  of  the  barrier  that  were  previously  combined  through  transmission  and  reflection 
at  the  barrier.  In  this  case, 


=  T(pj)f{xz+i,-g-j)  -  R(pj)fZ-j 
:+1J  T(pj )  -  R(Pj) 

=  T(pj)f{xz,-q-j)  -  R(Pj)fz+ i,-j 

Iz'j  T(pj)  -  R(pj ) 


for  j  >  0, 
for  j  <  0. 
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In  both  cases,  the  densities  f(xz+i,iq±j)  and  f(xz,±q±j)  are  approximated  in  a 
manner  similar  to  (4.11). 

To  approximate  to  second-order  in  the  p-direction  we  have 

1/2  =  "F  2  —  ^i)^Pap(fij) 

with  A i  —  \dxVi\At/ Ap  and  the  slope  erp(-)  defined  using  the  van  Leer  limiter  (4.12). 

5.  Numerical  examples.  In  this  section,  we  present  a  few  examples  of  both 
pure  state  dynamics  and  mixed  state  dynamics  in  order  to  verify  and  validate  the 
semiclassical  model  and  numerical  scheme. 

For  a  mixed  state  solution  with  a  macroscopic  distribution,  we  are  not  limited  by 
the  support  of  the  wavepacket  and  the  complexity  of  the  scheme  is  0((Aa:ApAt)_1) 
where  Ax,  A p,  and  At  3>  £■  For  direct  simulation  of  the  von  Neumann  equation,  not 
only  must  we  resolve  e  in  space  and  time  but  we  must  solve  the  equation  over  two 
space  dimensions  and  one  time  dimension  so  the  complexity  of  the  scheme  is  0(e-3). 
When  e«1,  the  computing  time  for  a  direct  von  Neumann  solution  is  considerably 
longer  than  for  the  multiscale  semiclassical  Liouville  solution. 

The  numerical  Schrodinger  solution  may  be  computed  using  the  Crank-Nicolson 
operator 


(5.1)  ip(xi,  t  +  At)  =  (1  +  ie  1AtHp)  1(l—ie  1  At  i/c)^  (*»>*) 


where  the  discrete  Hamiltonian  operator 


(5.2) 


Hd  = 


— e2  5j'j- 1  —  2<5jj  +  6jtj- )-i 

2m  (Ax)2 


+  V(Xi) 


with  Kronecker  delta  Su  =  1  and  5tj  =  0  if  i  ^  j.  Markowich  ef  al.  [24]  show  that  for 
such  a  scheme,  in  order  to  guarantee  correct  approximation  to  physical  observables 
for  small  e,  one  needs  to  take  Ax  =  o{e)  and  At  =  o(e).  One  may  also  compute 
the  numerical  Schrodinger  solution  using  a  pseudospectral  method  with  Strang  split¬ 
ting  [4].  In  this  case,  one  splits  the  kinetic  and  potential  terms,  so  that  for  each  time 
step 


■tp{x,  t  +  At)  =  eAtB/2T^  [eAM T  [eAtB/2^(x,  t)]  ] 

where 

A  =  —  k2  and  B  =  -  V(x) 

2  mi  ie 

and  the  operators  !F  and  denote  the  one-dimensional  discrete  Fourier  transform 
and  discrete  inverse  Fourier  transform  with  respect  to  the  x  and  k  variables.  One  can 
use  a  mesh  that  is  coarser  than  the  mesh  required  by  a  finite-difference  method  to 
resolve  e  and  capture  the  correct  dynamics  [4,  5].  Based  on  numerical  observation, 
we  find  that  we  require  Ax  <  e/4  to  ensure  numerical  convergence  to  the  correct 
physical  observables  and  for  numerical  error  to  be  insignificant.  When  the  potential 
is  discontinuous,  we  find  that  the  solution  exhibits  artificial  oscillations  unless  At  < 
(Ax)2/e  and  A t<e/V(x). 

The  von  Neumann  equation 


=  Hp-  ( HpT)T 


with  H  =  dxx  +  V(x) 
2m  w 
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has  the  formal  solution 


p(x,  x'  ,t  +  At)  =  etEAt^  p(x,  x',  t)e  ieAt^ . 

By  using  the  discrete  Hamiltonian  operator  (5.2),  we  may  approximate  the  von  Neu¬ 
mann  solution  in  terms  of  the  Crank-Nicolson  operator  (5.1)  to  get  a  method  without 
splitting  error 

Pi*1  =  (1  +is~1AtHD)~1(l  -  ie~A AtH o)p*ji  with 

p*j  =  (1  -  i£~l  AtHo)~l  (1  +  i£~1AtHD)p'ji 

where  py  =  p(xi,x'j,tn).  We  may  also  solve  the  von  Neumann  equation  using  a 
pseudospectral  method  with  Strang  splitting  [13], 

pn+\  =  gAtB/2^-1  f&tAj:  ^eMB/2pny^ 

where 

A  =  -~-.{k2  -  ka)  and  B  = -(V{x) -V(x')) 

2  mi  ie 

and  the  operators  T  and  T~x  denote  the  two-dimensional  discrete  Fourier  transform 
and  discrete  inverse  Fourier  transform  with  respect  to  the  (x,  x')  and  (k,  k')  variables. 
The  FFTs  may  be  optimized  by  exploiting  the  hermicity  of  the  density  matrix. 

Alternatively,  we  may  calculate  the  von  Neumann  solution  indirectly  by  solving 
the  Schrodinger  equation  for  several  states  and  then  using  definition  (2.7)  to  construct 
the  density  matrix.  This  simplifies  a  two-dimensional  problem  over  N2  grid-points  to 
n  independent  one-dimensional  problems  over  N  grid-points.  If  the  initial  distribution 
is  localized  in  phase  space,  n  may  be  chosen  to  be  appreciably  smaller  than  N,  saving 
not  only  memory  but  also  contributing  to  a  considerable  reduction  in  computation 
time.  Furthermore,  this  approach  allows  us  to  implement  the  solution  using  a  parallel 
computer  cluster.  One  way  to  implement  such  a  scheme  is  to  use  states  generated  by 
taking  thin  slices  of  the  initial  distribution  along  the  T-direction.  Consider  the  WKB 
initial  condition 

(5.3)  ip{x,Q\x,p)  =  {ax\/2T:)~1^  —  x)2 /4.a2)£xp{ipx/e), 

which  describes  a  wave  packet  with  an  0(1)  spread  in  position  and  O(e)  spread  in 
momentum.  Let  the  weight  distribution  in  the  definition  of  the  density  matrix  (2.7) 
be 


f(x,p)  =  8(x  -  x0)  exp(— (p  —  p0)2/2s2eal)/{s2apV2n) 


where  the  scaling  factor  s£  =  1/ i/l  +  (e/2 crxap)2.  Then 
p(x,  x1 ,0)  =  JJ  f(x,p)ijj(x,0-,x,p)ip(x',0;x,p)dxdp 


f  (x-  x0)2  +  ( x '  -  xq )2  (x  -  x')2  ipo{x  -  x') 

:eXP(  4ct|  2  s2s7W  £ 


1  /  {h(x  +  x>)  ~xo)2  ( x-x ')2  ip0(x-x') 

axV 2^eXPV  2a2  ~  2eV2  e 
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Using  the  Wigner  transform  (2.10),  we  have  the  equivalent  Liouville  initial  distribu¬ 
tion 


(5.5) 


f(x,p,  0)  = 


1  /  (x  —  Xp)2 

2 ir<jx(Tp  6XP  V  2 a2 


(p-Po)2\ 

H  J 


which  is  independent  of  e. 

To  compare  the  convergence  of  the  Schrodinger  and  von  Neumann  solution  in  the 
semiclassical  limit,  we  use  the  L1-error  of  the  position  probability  density  function 

(pdf). 


/OO 

\p(x,t)  -  \ip(x,  f)|2|  dx 

-OO 

with  p(x,t )  =  J_  f(x,p,t)  dp.  We  replace  \ip(x,t)\2  with  p(x,x,t)  for  the  von  Neu¬ 
mann  solution.  The  semiclassical  Liouville  model  should  also  predict  the  correct  weak 
limit  for  multiphase  solutions  when  interference  in  the  Schrodinger  and  von  Neumann 
solutions  result  in  oscillations  in  the  probability  density  distribution.  To  measure 
the  weak  convergence  in  the  semiclassical  limit,  we  determine  the  L1 -error  in  the 
cumulative  distribution  function  (cdf),  i.e.,  the  antiderivative  of  position  density  [11] 


p(s,t)  -  \ip{s,  f)|2  ds 


dx. 


In  each  example  we  compare  the  exact  or  numerical  semiclassical  Liouville  so¬ 
lution  with  numerical  Schrodinger  or  von  Neumann  solutions  for  equivalent  initial 
distributions  and  potentials.  Since  the  interactions  with  the  boundaries  are  not  rel¬ 
evant  to  the  study,  a  sufficiently  large  domain  is  chosen  and  simulation  is  stopped 
before  the  wave  envelope  reaches  the  boundaries. 

5.1.  Schrodinger  0(1)  wave  envelope  with  a  step  potential.  Consider  the 
step  potential 


(5.6) 


if  x  <  0, 
if  x  >  0. 


A  particle  impinging  on  this  potential  from  the  left  is  totally  reflected  when  the 
incident  velocity  is  less  than  1. 

We  find  the  exact  solution  by  the  method  of  characteristics  by  tracing  along  the 
bicharacteristics  backward  in  time  to  the  initial  conditions.  Let  0(t)  =  {  (x,p)  |  x  < 
0  and  x  —  pt  <  0,  or  a:  >  0  and  x  —  pt  >  0  }  be  the  region  in  phase  space  where  the 
bicharacteristics  have  not  crossed  the  quantum  barrier  at  x  =  0  within  a  time  t.  Then 
the  exact  solution 


(5.7)  f{x,p,t) 


f(x-pt,p,0),  (x,p)en{t) 

T  ■  f  x  —  qt,q,  0^  +  R  ■  f  (— x  +pt,  -p,  0) ,  otherwise 


where  the  incident  velocity  is  given  by  q  =  \/p2  + 1  if  p  >  0  and  q  =  —yjp2  —  1  if 
p  <  0.  From  equations  (4.3)  and  (4.7),  the  reflection  coefficient  is  given  by 


R  = 


p-q 

p  +  q 


\p  ~  <?l4- 
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Note  that  when  p  €  [—1, 0],  q  is  imaginary  and  R  =  1  indicating  total  reflection. 
Consider  the  WKB  initial  condition 

il>(x,0)  =  A(x)eiS^E 

as  a  wavepacket  generalization  with  the  amplitude  and  phase  functions  given  by 

A(x)  =  (7T  fT2)-l/4e-(3:-xo)2/CT2 
S  (x)  —  ax 2  +  bx  +  c. 

Since  S(x)  is  a  quadratic  function,  we  can  calculate  the  Wigner  transform  of  ip(x,  0) 
exactly  to  get 


f(x,p,  0)  =  (■Ke)-1e(-x-Xo)2/a2e-^  ax+b-pf/^/^)\ 

In  the  semiclassical  limit  (e  — >  0),  we  have 

f(x,p,  0)  =  A2(x)S(p  -  'VxS(x)) 

(5.8)  =  (ay/Tt)~1e~(x~x°^2/a26(p  -  (2 ax  +  b)). 

By  taking  a  =  0(1)  in  A(x),  we  create  a  wave  envelope  that  is  independent  of 
£,  allowing  us  to  study  the  convergence  of  solutions  as  e  — ►  0.  When  a  ^  0,  the 
distribution  of  phases  included  in  the  Schrodinger  solution  is  also  0(1). 

Using  the  above  semiclassical  WKB  initial  conditions  (5.8),  we  note  that  when 
t  =  —  l/2a  the  position  density  for  the  Liouville  solution  (5.7)  exhibits  a  caustic 
with  all  bicharacteristics  intersecting  at  either  x  =  b/2a,  or  x  —  —b/2a  for  reflected 
solutions.  Because  of  the  nonlinear  change  to  the  incident  velocities,  the  transmitted 
bicharacteristics  do  not  cross  simultaneously  resulting  in  a  traveling  front,  the  leading 
edge  of  which  is  unbounded. 

Take  (zo,po)  =  (—5, 1)  and  take  a  =  —  b  =  po  —  2a  =  |  and  a  =  -^j.  Then  we 
have  the  initial  conditions 


ip(x,  0)  =  (10/7r)1/4e-2°o(a:-!Co)2ei(a:c2+(Po-2o*o)®)6: 


for  the  Schrodinger  equation  and 

f{x,p,  0)  =  (10/7T)1/2  e-100(T-Io)2j(p  _  po  -  2a(x  -  x0 )) 

for  the  semiclassical  Liouville  equation.  The  numerical  Schrodinger  solution  is  solved 
using  a  Crank-Nicolson  finite-difference  method  over  the  domain  [—1,1]  with  mesh 
size  Ax  =  A t  =  10~6.  The  exact  semiclassical  Liouville  solution  is  solved  by  tracking 
characteristics  forward  in  time  with  values  determined  by  the  initial  velocity  given  by 
VS  =  |  —  5 x .  We  compute  the  solution  at  time  t  —  0.8. 

The  position  densities  for  several  values  of  e  are  shown  in  Fig.  5.1.  The  conver¬ 
gence  results  of  the  errors  in  the  two  solutions  are  listed  in  Table  5.1.  Based  on  this 
study,  we  find  that  the  l 1  convergence  rate  in  e  of  the  pdf  is  about  0.6  and  the  l1 
convergence  rate  in  e  of  the  cdf  is  about  1.1. 
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Table  5.1 

Errors  in  solutions  for  different  values  of  £  for  Example  5.1. 


£ 

200 -1 

800“ 1 

3200- 1 

12800- 1 

i1 -error  (pdf) 

8.78  X  10-1 

3.37  x  10"1 

1.55  x  10-1 

8.61  x  lO-2 

l1  -error  (cdf) 

5.15  X  10~2 

1.00  x  10“2 

2.28  X  10-3 

1.08  x  10~4 

FlG.  5.1.  Position  densities  for  the  semiclassical  Liouville  (top)  and  Schrodinger  (bottom) 
solutions  of  Example  5.1.  The  Schrodinger  solution  shows  e  =  (a)  200“ 1,  (b)  800-1,  (c)  3200-1 
and  (d)  12800-1.  The  position  density  of  Liouville  solution  exhibits  a  caustic  near  x  =  0.08  and  the 
peak  is  unbounded.  For  the  Schrodinger  solution  the  peak  reaches  a  height  of  19  for  the  £  =  12800-1 . 
The  plots  are  truncated  for  clarity. 


5.2.  Von  Neumann  solution  with  step  potential.  We  now  consider  the 
solution  to  the  von  Neumann  equation  with  the  step  potential  given  Example  5.1. 
To  construct  a  von  Neumann  initial  condition  p{ x,x',0)  which  corresponds  to  a  Li¬ 
ouville  initial  condition  f(x,p,  0),  we  may  directly  use  the  definition  of  the  density 
matrix  (2.7)  for  some  weight  function  with  the  probability  amplitudes  tp(x,t)  given 
by  Gaussian  e-wavepackets 

(5.9)  l/.(z,  0)  =  (7r£)-l/4e-(x-xo)2/2eeipox/e_ 

The  Liouville  initial  condition  may  subsequently  be  calculated  by  a  Wigner  transform 
of  the  density  matrix.  Alternatively,  we  may  construct  the  density  matrix  by  using 
the  inverse  Wigner  transform  applied  to  the  Liouville  initial  conditions  f(x,p,0)  to 
get 


/OO 

f(Hx  +  x')’Pi  0)e!p(l_I')/£  dp. 

-OO 
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By  taking  the  Liouville  initial  conditions  to  be  the  Gaussian 


(5.10) 


f{x,p,Q)  = 


1 


■  exp 


— (s  -  xq): 


)  exp  ^ 


-ip-Pof 


K 


2naxap  ~ V  2a2 

we  may  compute  the  von  Neumann  initial  conditions  exactly  to  get 


(5.11)  p(x,x',Q)  = 


1 


ox\f2n 


exp 


(5(2:  +  x')  -  rc0)2  {x  -  x')2  ip0(x  —  x') 


2  ol 


2 e2av  2 


We  chose  ax  =  ap  =  0.05,  xq  =  —0.5  and  po  =  1.0  and  compared  the  solutions  to 
the  von  Neumann  and  semiclassical  Liouville  equations  at  time  t  =  1.0.  The  von  Neu¬ 
mann  equation  was  solved  using  the  psuedospectral  method  with  Strang  splitting  over 
the  domain  [—1,1]  with  e  =  64_1,  128-1,  256-1  and  512-1.  The  grid  spacing  was 
fixed  at  Aa;  =  2048 _1  with  At  —  [Ax)2  /e  to  ensure  consistency  and  stability.  The 
exact  solution  to  the  semiclassical  Liouville  model  is  calculated  using  equation  (5.7). 

The  position  densities  for  the  semiclassical  Liouville  solution  and  the  von  Neu¬ 
mann  solution  for  several  values  of  e  are  shown  in  Fig.  5.2.  The  errors  in  the  two 
solutions  are  listed  in  Table  5.2.  Based  on  our  study,  we  find  the  convergence  rate  of 
the  ^-error  of  the  pdf  is  about  0.7  as  e  —>,0  and  the  convergence  of  the  /1-error  of  the 
cdf  is  about  0.9  as  e  — >  0. 


Table  5.2 

Errors  in  solutions  of  Example  5.2  for  different  values  of  e  . 


e 

64“ 1 

128-1 

256- 1 

512“ 1 

P-error  (pdf) 

6.03  x  10-1 

4.04  x  10-1 

2.50  x  Kr1 

1.40  x  10-] 

l1 -error  (cdf) 

9.22  X  10~2 

4.83  x  10-2 

2.53  x  lO-2 

1.32  x  lO"2 

We  may  also  consider  the  effect  of  incorporating  barrier  time  delay  in  the  ap¬ 
proximation  of  the  von  Neumann  equation  for  nonvanishing  e.  As  evident  from  the 
offset  of  the  centers  of  the  distributions  on  the  left  side  of  Fig.  5.2,  one  source  of 
error  is  the  time  delay  which  vanishes  in  the  semiclassical  limit.  The  time  delay  may 
be  considered  as  an  0(e)  correction  and  hence  we  may  neglect  it  in  the  semiclassical 
limit.  While  the  the  addition  of  a  delay  time  is  numerically  nontrivial,  for  the  analytic 
solution  (5.7)  it  is  a  straight-forward  modification. 

Typically,  time  delay  is  defined  in  terms  of  Wigner  time  delay,  the  delay  to  the 
group  velocity  of  a  wave  packet  resulting  from  reflection  and  transmission.  As  such 
it  is  meaningful  when  the  wave  packet  has  a  well-defined  peak.  This  is  not  generally 
the  case,  especially  when  the  barrier  is  sufficiently  wide.  Considering  the  scattering 
relation  (4.5),  the  reflection  and  transmission  group  delay  times  for  unit  mass  are  [26] 


e  d  eT  .  1  dt 

n  =  - ~r  arg t  =  -Im(-— ) 
p  dp  p  t  dp 


e  d  e  1  dr 

and  Tr  =  argr  =  -Im  (-— ). 
p  dp  p 


V  dp 


For  the  step  potential  (5.6),  we  have  from  equation  (4.3)  that  the  reflection  time 
delay  is 


t>  =  2elm  [(pq)  x]  =  ~  7= — = 

pyi  ~P 


when  p  €  [—1,0].  There  is  no  transmission  or  reflection  delay  time  for  p  £  [—1,0].  To 
incorporate  the  time  delay,  we  make  the  replacement 


f{-x+pt,-p,  0)  ->  f  (-X+p(t  +  Tr),-p,  0) 
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X 


Fig.  5.2.  Position  densities  for  the  semiclassical  Liouville  (top)  and  von  Neumann  (bottom) 
solutions  of  Example  5.2.  The  von  Neumann  plot  shows  e  equal  to  (a)  64_1,  (b)  128-1,  (c)  256_1 
and  (d)  512"1. 


in  the  reflected  term  of  the  exact  solution  (5.7). 

We  compare  the  von  Neumann  solution  with  the  Liouville  solution  with  time  delay 
correction.  The  ^-errors  are  listed  in  Table  5.3.  Based  on  this  study,  we  find  that  the 
addition  of  delay  time  provides  some  improvement  to  the  model.  The  convergence 
rate  of  the  Z1-error  of  the  pdf  is  about  1.3  and  convergence  rate  of  the  ^-error  of  the 
cdf  is  about  0.9  as  e  — >  0. 

Table  5.3 

Errors  in  solutions  of  Example  5.2  with  time  delay  correction. 


e 

64-1 

128“ 1 

1 

CO 

in 

512-1 

-error  (pdf) 

3.67  x  io-1 

1.78  x  10-1 

7.05  X  IO"2 

2.23  x  IO"2 

i1 -error  (cdf) 

2.62  x  10-2 

1.65  x  10“2 

7.80  x  IO"3 

3.90  x  10-3 

5.3.  Von  Neumann  solution  with  two  step  potentials.  We  may  consider 
more  complicated  geometries  by  considering  multiple  barriers.  In  this  example  we 
construct  an  0(l)-wide  rectangular  barrier  by  taking  two  step  barriers  sequentially. 
Consider 


if  Z  €  [0,  |], 
otherwise. 


We  take  the  initial  conditions  given  in  equations  (5.10)  and  (5.11)  with  crx  =  crp  =  0.05, 
xq  —  —0.45  and  po  =  1.1.  We  compute  over  the  domain  [—1.25, 1.25]  and  compare  the 
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solutions  at  time  t  =  1.2.  The  von  Neumann  equation  is  solved  using  a  pseudospectral 
method  with  Strang  splitting  as  in  Example  5.2.  The  semiclassical  Liouville  solution 
is  solved  using  the  numerical  method  proposed  in  §4  using  N  grid  points  in  x  and 
p  and  1.57V  steps  in  time.  The  results  are  shown  in  Fig.  5.3  with  e  =  0.002.  Even 
with  a  fairly  coarse  mesh,  the  numerical  semiclassical  solution  agrees  well  with  the 
von  Neumann  equation  both  in  the  strong  limit  away  from  the  barrier  and  in  the  weak 
limit  between  the  two  step  potentials.  See  Fig.  5.4. 

We  calculate  convergence  rate  as  Ax ,  A p,  At  — >  0  of  numerical  scheme  for  the 
semiclassical  Liouville  equation  by  computing  the  terror  of  the  numerical  solutions 
using  a  mesh  with  N  =  50,  100,  200,  and  400  grid  points.  For  an  “exact”  solution, 
we  use  the  numerical  solution  using  N  =  3200.  The  errors  are  listed  in  Table  5.4. 
Based  on  this  study,  we  find  the  convergence  rate  of  the  numerical  scheme  using  the 
Z1-norm  is  about  1.2. 

Table  5.4 

Errors  in  solutions  of  Example  5.3  for  various  mesh  sizes  Ax. 


grid  points 

50 

100 

200 

400 

l1  -error 

3.32  x  10-1 

1.15  x  10-1 

4.72  x  10~2 

2.56  x  lO-2 

5.4.  Resonant  tunneling  von  Neumann  solution.  We  present  a  final  exam¬ 
ple  to  illustrate  a  specific  physical  model,  the  resonant  tunneling  diode  (RTD)  [19, 
25,  30],  An  RTD  consists  of  thin  layers  (a  few  nanometers  thick)  of  different  semicon¬ 
ductors,  such  as  gallium  arsenide  (GaAs)  and  aluminum  gallium  arsenide  (AlGaAs), 
that  are  sandwiched  together  to  form  a  double-barrier  quantum  well  structure.  For 
semiconductors  the  de  Broglie  wavelength  is  on  the  order  of  tens  of  nanometers,  so  the 
length  of  the  entire  RTD  structure  is  on  the  length  scale  of  a  de  Broglie  wavelength. 
The  region  outside  the  barrier  is  doped  to  provide  a  sufficient  number  of  free  electrons. 
Unlike  the  transmission  probabilities  of  the  step  potentials  presented  the  previous  ex¬ 
amples,  the  transmission  probability  of  an  RTD  is  not  a  monotonic  function  of  the 
incident  particle  energy.  Rather,  it  is  oscillatory  and  admits  narrow  peaks  of  total  or 
almost  total  transmission  well  below  the  cutoff  energy  for  classical  transmission.  By 
changing  the  bias  voltage  of  an  external  electrostatic  potential  applied  to  the  system, 
the  resonance  may  be  tuned  to  admit  electrons  of  varying  energies. 

We  shall  assume  that  the  electron  trajectory  is  ballistic.  In  the  quantum  region, 
this  simplification  is  appropriate  since  the  electron  mean  free  path  is  substantially 
larger  than  the  barrier  thickness.  However,  away  from  the  barrier  this  simplification 
is  physically  unrealistic  since  the  electron  mean  free  path  is  small  compared  to  the 
classical  length  scale  for  a  dense  medium.  In  this  case,  a  relaxation  term  or  collision 
operator  should  be  added  to  the  Liouville  equation  to  capture  the  particle  dynamics. 
Since  we  require  that  the  Hamiltonian  be  only  locally  preserved,  the  model  may  be 
extended  to  a  dissipative  system,  for  which  the  Hamiltonian  is  continuous,  without 
changing  the  approach  discussed  in  §3  and  §4.  Hence,  for  the  purpose  of  validation, 
the  assumption  is  reasonable. 

€  (— oo,  —a  —  b ] 

€  (—a  —  b ,  —  a]  U  (a,  a  4-  fe] 

€:  [ — a,  a  +  6] 

G  (a  +  6,  oo) 


We  construct  a  representative  barrier 

r+fvb 

5  Vo  x/(a  +  b)  +  V , 
\Vqx/  (a  +  b) 

i-IVo 


V(x)  =  { 
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Fig.  5.3.  Position  densities  for  the  numerical  semiclassical  Liouville  (top)  and  von  Neumann 
(bottom)  solutions  of  Example  5.3.  The  •  in  the  Liouville  plot  shows  the  numerical  solution  for  with 
150  grid  points  over  the  domain  [—1.25, 1.25].  The  solid  line  shows  the  numerical  solution  for  3200 
grid  points.  The  von  Neumann  solution  is  for  e  =  0.002. 


FlG.  5.4.  Detail  of  Fig.  5.3  showing  position  densities  for  the  numerical  semiclassical  Liouville 
and  von  Neumann  solutions.  The  •  shows  the  numerical  solution  for  with  150  grid  points  over  the 
domain  [—1.25,1.25].  The  solid  line  shows  the  “exact”  Liouville  solution  and  the  von  Neumann 
solution  using  e  =  0.002. 
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FlG.  5.5.  Transmission  probability  as  function  of  the  momentum  p  for  the  RTD  barriei — shown 
in  the  inset — presented  in  Example  5.4 . 


where  the  external  potential  bias  Vq  =  0.48,  the  thickness  of  each  barrier  b  =  O.Qe,  the 
thickness  of  the  well  separating  the  barriers  2 a  =  1.2s,  and  the  height  of  each  barrier 
V ,  —  2.25.  See  Fig.  5.5.  We  take  Gaussian  initial  distributions  (5.10)  and  (5.11)  with 
ax  =  0.05,  ap  =  0.15,  xq  =  -1  and  p0  =  1.  The  solutions  are  computed  over  the 
domain  [—4, 4]  and  compared  at  time  t  =  2.5. 

The  von  Neumann  equation  is  solved  indirectly  using  the  WKB  initial  condi¬ 
tions  (5.3)  with  weight  distribution  (5.4).  We  use  a  Crank-Nicolson  finite-difference 
method  to  solve  the  Schrodinger  equations.  To  ensure  that  the  weight  function  is 
sufficiently  resolved,  we  take  N  =  (5e)-1  Schrodinger  solutions  with  initial  values 
equally  spaced  over  8ap  about  po. 

The  semiclassical  Liouville  solution  is  solved  using  the  numerical  method  proposed 
in  §4  using  an  N  grid  points  over  [—4, 4]  in  x,  2 N  grid  points  over  [—3, 3]  in  p  and  3 N 
steps  in  time.  The  exact  solution  is  computed  using  equation  (5.7)  with  transmission 
and  reflection  probabilities  calculated  using  the  transfer  matrix  method.  In  computing 
the  transfer  matrix  for  both  the  numerical  and  exact  solutions,  the  quantum  barrier 
is  discretized  using  2000  grid  point  over  the  length  6e  for  an  arbitrary  e.  The  results 
are  shown  in  Fig.  5.6  and  Table  5.5.  We  calculate  an  ^-convergence  rate  of  1.7  in 
Ax,  Ap,  At. 


Table  5.5 

Errors  in  solutions  of  Example  5.4  for  various  mesh  sizes  Ax. 


grid  points 

80 

160 

320 

640 

l1 -error 

3.01  X  10-1 

1.37  x  10"1 

4.43  x  10~2 

8.90  x  10~3 
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FlG.  5.6.  Position  densities  for  the  numerical  semiclassical  Liouville  (top)  and  von  Neumann 
(bottom)  solutions  of  Example  5.4 •  The  Liouville  solution  shows  the  numerical  solution  for  (a)  80, 
(b)  160, (c)  320  and  (d)  640  grid  points.  The  von  Neumann  solution  is  plotted  for  e  =  (a)  50-1, 
(b)  100”1 11,  (c)  200-i  and  (d)  400-1.  The  exact  solution  (*)  is  also  plotted  in  each  case. 
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